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Abstract 


The framework for constructing a high-order, conservative Spectral (Finite) Volume (SV) 
method is presented for two-dimensional scalar hyperbolic conservation laws on unstructured 
triangular grids. Each triangular grid cell forms a spectral volume (SV), and the SV is further 
subdivided into polygonal control volumes (CVs) to supported high-order data reconstructions. 
Cell-averaged solutions from these CVs- are used to reconstruct a high order polynomial 
approximation in the SV. Each CV is then updated independently with a Godunov-type finite 
volume method and a high-order Runge-Kutta time integration scheme. A universal 
reconstruction is obtained by partitioning all SV. in a geometrically similar manner. The 
convergence of the SV method is shown to depend on how a SV is partitioned. A criterion based 
on the Lebesgue constant has been developed and used successfully to determine the quality of 
various partitions. Symmetric, stable, and convergent linear, quadratic, and cubic SVs have been 
obtained, and many different types of partitions have been evaluated. The SV method is tested for 
both linear and non-linear model problems with and without discontinuities. 

Key Words: Fligh-Order, Unstructured Grid, Spectral Volume, 2D Conservation Laws 
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1. Introduction 


We continue the development of the spectral (finite) volume (SV) method for hyperbolic 
conservation laws on unstructured grids following the one-dimensional framework presented in 
[38], We wish to pursue a numerical method for conservation laws which has all of the following 
properties: a) conservative, b) high-order accuracy, i.e., the order of accuracy is greater than 
second order, c) geometrically flexible, i.e., applicable for unstructured grids, and d) 
computationally efficient. The SV method is developed to hopefully satisfy these four 
requirements, in a relative sense with respect to the current state-of-the-art numerical methods 
such as the high-order k-exact finite volume (FV) method [5,16], essentially non-oscillatory 
(ENO) method [1,10,19], and weighted ENO (WENO) method [3,17,20,22,27], and the 
discontinuous Galerkin (DG) method [2,6,12-14], amongst many others. 

One of the most successful algorithms for conservation laws is the Godunov method [18], which 
laid a solid foundation for the development of modem upwind schemes including MUSCL [36], 
PPM [15], ENO [19] and WENO schemes [22,27], There are two key components in a Godunov- 
type method. One is data reconstruction, and the other is the Riemann solver. The original 
Godunov scheme employed a piece-wise constant data reconstruction, and the exact Riemann 
solver, and the resultant scheme was only first-order accurate. Later on, higher-order polynomial 
reconstructions are used to replace the piece-wise constant reconstruction, and approximate and 
more efficient Riemann solvers [21,25,29,31,35,37] were employed to substitute the exact 
Riemann solver. In addition, limiters were also introduced to remove spurious numerical 
oscillations near steep gradients [36] in higher than first-order Godunov method. 

Although Godunov-type methods were originally developed for structured grids, they have been 
successfully extended to unstructured grids, thus achieving greater geometric flexibility. Most of 
the unstructured grid methods are second-order accurate because they are relatively easy to 
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implement, and are quite memory efficient. Several high-order schemes have been developed for 
unstructured grids. For example, a high-order k-exact finite volume scheme was developed by 
Barth and Fredenckson in [5], An ENO scheme for unstructured grid was developed by Abgrall 
in [1]. Two WENO schemes for unstructured grids were developed by Fnednch in [17], and Hu 
and Shu in [20], Although high-order accurate finite volume schemes can be obtained 
theoretically for an unstructured grid by using high-order polynomial data reconstructions, higher 
than linear reconstructions are rarely used in three dimensions in practice. This is mainly because 
of the difficulty in finding valid high-order (non-singular) stencils, and the enormous memory 
required to store the coefficients used in the reconstruction. In a k-exact finite volume method, 
each control volume has a different reconstruction stencil. Therefore, a data reconstruction must 
be performed at each iteration for each control volume. This reconstruction step is the most 
memory and time consuming in higher than second-order schemes. In a recent implementation of 
a third-order FV scheme with a quadratic reconstruction in three dimensions by Delanaye and 
Liu [16], the average size of the reconstruction stencils is about 50-70. Still there are many 
singular reconstruction stencils. The size of the reconstruction stencils usually increases non- 
linearly with the order of accuracy. For a fourth order FV scheme, the average stencil size is 
estimated to be at least 120. It is very memory and CPU intensive to perform the reconstruction. 

More recently, another high-order conservative algorithm called the Discontinuous Galerkin 
method was developed by Cockbum, Shu, et al in a senes of papers [12-14] and also [2, 6] on 
unstructured gnds. In the DC method, a high-order data distribution is assumed for each element. 
As a result, the state variable is usually not continuous across element boundanes. The fluxes 
through the element boundanes are computed using an approximate Riemann solver, similar to 
FV methods. The residual is then minimized with a Galerkin approach. Due to the use of 
Riemann fluxes cross element boundanes, the DG method is fully conservative. A disadvantage 
of the DG method is that high-order surface and volume integrals are necessary, which can be 
expensive to compute. Another high-order conservative scheme for unstructured quadnlateral 
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grids is the multi-domain spectral method on a staggered grid developed by Kopriva et al [23- 
24], The multi-domain spectral method is similar to the spectral element method by Patera [30], 
which is not conservative. Other related methods include cell-average based spectral method [9] 
and spectral element type method [34]. Although very high-order of accuracy was achievable 
with these methods, the methods are difficult to extend to other cell types such as triangles, or 

tetrahedral cells. 

In [38], the first paper in the series, a new conservative high-order SV method is developed for 
conservation laws on one-dimensional unstructured grids. Through the use of spectral volumes 
(SVs) and control volumes (CVs), the method is not only conservative, but very efficient as well. 
In this paper, we extend the SV method to two dimensions. In the next section, we first review 
the basic framework of the SV method on triangular gnds. In addition, we present a TVD Runge- 
Kutta time integration scheme. In Section 3, the reconstruction problem based on CV-averaged 
solutions is studied, and it is shown that the reconstruction problems on all triangles with a 
similar partition are identical. In Section 4, convergent reconstructions for high-order SV 
schemes are discussed, and the partition of a SV is shown to affect the convergence of the 
method. Section 5 discusses issues related to discontinuity-capturing and several TVD and TVB 
limiters are presented. In Section 6, numerical implementations of the SV method for both linear 
and non-linear scalar conservation laws are carried out, and accuracy studies are performed for 
both linear and non-linear wave equations to verify the numerical order of accuracy. The shock- 
capturing capability of the method is also demonstrated with the Burger's equation. Finally, 
conclusions and recommendations for further investigations are summarized in Section 7. 

2. Review of the Spectral Volume Method 

Consider the following multi-dimensional scalar conservation law: 
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on domain Qx[OJ] and Q c R 2 with the following initial condition 


(2.1b) 


w(x, yfi) — y ) ’ 

and appropriate boundary conditions on dQ . In (2.1), x and y are the Cartesian coordinates and 
(jt,y)eQ, re [0,7] denotes time, u is a state variable and / and g are fluxes in * and y 
directions, respectively. Domain Q is discretized into / non-overlapping triangular cells. In a k- 
exact FV method, a data reconstruction is performed for each cell using data from a collection of 
neighboring cells, collectively known as a reconstruction stencil as shown in Figure la. A unique 
flux through each face is then computed given the reconstructed state variables at both sides of 
the face using either an exact or approximate Riemann solver. A summation of fluxes through all 
the faces of a cell is then used to update the cell-averaged state variable. 

In the SV method, the triangular cells are called spectral volumes, denoted by S„ which are 
further partitioned into subcells named comrol volumes (CVs), denoted by C u . Volume-averaged 
state variables on the CVs are used to reconstruct a high-order polynomial inside the SV. To 
represent the solution as a polynomial of degree m in 2D, we need N = (m+l)<m+2)/ 2 pieces of 
independent information, or degrees of freedom (DOFs). The DOFs in a SV method are the 
volume-averaged mean vanables at the N CVs. For example, a SV supporting a quadratic data 
reconstruction is shown in Figure lb. Other candidate partitions for linear to cubic SVs are 
shown in Figures 2-4. The number of CVs in Figures 2-4 is the minimum required for these 
polynomial reconstructions. Other CV subdivisions are definitely posstble. Integrating (2.1) on 

Cij, we obtain 


du. 


(2.2) 


where F = (f g). and n is the unit outward normal of 8C U . the boundary of C,j. Define the CV- 
averaged state variable for C t j as 

| udV 


- -fk 

u uj y 


(2.3) 


ij 


where V u is the volume (area in 2D) of Cq. Then (2.2) becomes 


V i,j r= 1 A r 


(2.4) 


where K is the total number of faces in Qj, and A r represents the r-th face of C u . The surface 
integration on each face can be performed with a k-th order accurate Gauss quadrature formula 

(, k=m+l ) , i e. 


f (F • n)dA = w rq F(u(x rq y rq )) • n r A r + 0{A r h k ), 

i 


(2.5) 


where J = integer[(k+l)/2] is the number of quadrature points on the r-th face and, w rq are the 
Gauss quadrature weights, (x rq , y rq ) are the Gauss quadrature points, h is the maximum edge 
length of all the CVs. Time t is omitted whenever there is no confusion. If F = constant , the 

following identity exists: 


^ J(F • n)dA = 0. 


( 2 . 6 ) 


r= 1 A. 


Therefore, we will gain an extra order of accuracy if we sum up the surface integrals for the 
faces of C u , i.e.. 


V \(F*n)dA = j^f J w rq F(u(x rq 'y rq ))*n r A r +0(A r h k+ '). 

r=l A '• =1 " =1 


(2.7) 


Since 0(V,) = 0( A r h) , we therefore have 


Ly l f • n)dA = -k)) # nA + 0{hk y 

V U £ll V ur-l,-l 


( 2 . 8 ) 


6 


Now assume a multi-dimensional polynomial in * and y of order at most k - 1 exists on 5, which 
is a k-th order approximation to the state variable, i.e., 

p.(x,y) = u(x,y) + 0(h k ), (. x,y)e S r ( 2 - 9 ) 

With the polynomial distribution on each SV, the state variable is most likely discontinuous 
across the SV boundaries, unless the state variable is a polynomial of order k - 1 or less. 
Therefore, the flux integration involves two discontinuous state variables just to the left and right 
of a face of the SV boundary. This flux integration is carried out using an exact Riemann solver 
or one of the Lipschitz continuous approximate Riemann solvers or flux splitting procedures, i.e., 
F{u(x rq , y rq )) • n r = F Riem {Pi (x rq , y rq lPiA X rq ’ ) + 0(Pi(x rq ,y rq ) ~ PiA X rq ’ < 2 ' 10) 

Here Pi r is the reconstruction polynomial of a neighboring CV C u , r , which shares the face A r 
with C/j. Both Pi and p i r are k-th order approximations of the exact state variable, i.e.. 

Pi (-v>%> = u ^yA + °^ k )’ (2J la) 

PiAx rr y rq )-^x rr y rq )AO{h k ). (2Ub) 

Therefore 

F{u(x rq ,y rq ))*n r = F Riem (p i (x rq ,y rq ),p i Ax r<l ,y rq ),n r ) + 0 { hk ). < 2 ' 12) 

Substituting (2.12) into (2.4), we obtain 

f (F • n)dA = 'Zw rq F Rie Jp l (x r(l ,y rq ), p u r (x rq , y rq \ n , K + °^ k )• (2 ' 13) 

i 

Summarizing (2.4)-(2.13), we obtain the following semi-discrete, k-th order accurate scheme on 
Cij for the conservation law (2.1) 


du. 


'FJ 


- , =0(h ‘ x ai4) 

Clt V ijr=lq = 1 

For time integration, we will use the third-order TVD Runge-Kutta scheme from [32], We first 
rewrite (2.14) in a concise ODE form 




(2.15a) 
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where 



M l,l 


'Ru 

(w) 

ll - 

Kj 

^(«) = 


•(“) 


W /,;V _ 




and 


R u = — ^Yj W rq F R^ Pi ^ X r^ y P '^ Xrq,y r ^ ,tlr>)Ar ' 

ViJ r=l 9=1 

Then the third-order TVD Runge-Kutta scheme can be expressed as: 


« 0) =«" +&R h (u n )-, 


(2.15b) 


(2.15c) 


u {2) =-u n +-[m ( 1) +A/i? fc (u (1> )]; (2 ‘ 16) 

4 4 

U n+ 1 = -«" +-[W (2) + Arif*(u (2) )]. 

3 3 

The SV method idea can of course be easily extended to other cell types such as quadrilaterals, 
tetrahedra, hexahedra, prisms, etc. For cell types other than triangles and tetrahedra, it seems 
symmetric CV subdivisions with the minimum number of CVs for a given order of accuracy are 
difficult to obtain. The development of the SV method for other cell types will be reported 

elsewhere. 


3. The Reconstruction Problem in A Spectral Volume 


As discussed in the previous section, in order to compute the flux across a surface, one needs to 
evaluate the state vanable u at quadrature points. These evaluations can be achieved by 
reconstructing the state vanable u in terms of some basis functions using the CV-averaged 
solutions U j within a SV. (For simplicity, we drop the subscript i and use u tj , Cj and V; to 
denote C,, and V UJ , respectively.) In general, one can choose any linearly independent 
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functions as the basis functions. Here we focus only on the reconstruction using polynomials as 
the basis functions. 


N_ = 


Let P m denote the space of degrees/ polynomials in two dimensions. Then the dimension of the 
approximation space is 

m + 2^| ( m + 1 )(m + 2) 

2 = 2 ' 

/ 

which is the minimum dimension of the space that allows P„ to be complete. In order to 
reconstruct u in P„ we need to partition the SV into a set of N m nonoverlappmg CVs. Let S 
denote the physical space of a SV and n„, denote the partition, i.e., 

n m = {C V C 2 ,-,C W } , 

where C ■ C S,j = (or N if there is no confusion) is the j-th CV inside the 5 V, and 


s = U c ,- 

> = l 

The reconstruction problem reads: Given a continuous function in S, ueC{S) (the space of 
continuous functions in S), and a partition U m of 5, find p m e P m , such that 

f P m (x,y)dV =| u(x,y)dv, j = N m . (31) 

JCj j 

To actually solve the reconstruction problem, we introduce the complete polynomial basis, 
e, (x,y)eP m . where P. = spank, U y ) IS • Therefore p m can be expressed as 

p, = f>,(.t,y), < 3 - 2a > 


/ = 1 


or in the matrix form 

(3.2b) 

p m =ea, 

where e is the basis function vector 1 and a is the reconstruction coefficient vector 

[a, ,...,a N ] T ■ Substituting (3.2a) into (3.1), we then obtain 

‘M,yW=u,, j = (33) 


/=! 


Let if denote the column vector [if,, ] r , Eq. (3.3) can be rewritten in the matrix form 
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Ra = u, 

where the reconstruction matrix 


R = 


(3.4) 


[-Lj 

e, (x, y) dV • 

• fJ 

e N {x,y)dV 

V 

c, 


C ,1 

l 

f e l (x,y)dV • 


[ e N (x, y) dV 

V * 

v N 

lc jV 

V N ‘ 



(3.5) 


(3.6) 


The reconstruction coefficients a can be solved as 

a = R ] u, 

provided that the reconstruction matrix R is nonsingular. Substituting Eq. (3.6) into Eq. (3.2), 
p m is then expressed in terms of cardinal basis functions or shape functions L = [L, L v ] 

p^L^y^Lu. Q-7) 

i = i 


Here L is defined as 

L = eR~\ (3,8) 

Equation (3.7) gives the functional representation of the state variable u within the SV. The 
function value of u at a quadrature point or any point ( x rq , y rq ) within the SV is thus simply 

<3 ' 9) 

The above equation can be viewed as an interpolation of a function value at a point using a set of 
cell averaged values with each weight equal to the corresponding cardinal basis functional value 

evaluated at that point. 


Note that once the polynomial basis functions e, are chosen, the cardinal basis functions L j are 
solely determined by the partition U m of 5. The shape and the partition of 5, in general, can be 
arbitrary as long as the reconstruction matrix R is nonsingular. However, different shapes of 
spectral volumes can result in the same expression of the cardinal basis functions (in terms of a 
few geometric parameters) if a geometrically similar partition can be applied to them. In the 
following, we shall examine a special case that all SVs are triangular and all CVs are polygons 
with straight edges. In this case, even though the shapes of the SVs may all be different, as long 
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as they are partitioned in a geometrically similar manner, they all have the same reconstruction, 
in which the functional values of the cardinal bases at similar grid points are all exactly the same. 
We shall defer the discussion of other type of spectral volumes (e.g., quadrilateral, tetrahedron, 

curved boundaries, etc.) elsewhere. 


We first consider a transformation ¥ : S -> D , shown in Figure 5a, which transforms an 
arbitrary triangle 5 to a right triangle D. Another often-used transformation is from an arbitrary 
triangle to an equilateral triangle E, shown in Figure 5b. Let us use (x,y) to denote the 
coordinates in 5 and g.ij) the coordinates in D. For simplicity, we assume one of the nodes is 

located at the origin r 0 = (0,0) and the other two at r, =(*„?,) and r 2 =(x 2 ,y 2 ) m S, 
corresponding to (0,0), (1,0) and (0,1) in D, respectively. Thus, the transformation can be written 


y . r = r,£ + r 2 rj, 


£ >0,?7>0, 

and£ + T] < 1. 


(3.10) 


Since the transformation is linear, for a complete set of basis functions e(x,y)e P m , one can 


easily show that 

e(x,y) = e($,r])T. (3 * U) 

Here T is the transformation matrix containing only the geometric information of the nodal 
positions of 5. For example, if e(x,y) = [l,x,y,x 2 ,xy,y ], then 

‘1 0 0 0 0 0 

0 x ] y, 0 0 0 

0 x 2 y 2 0 0 0 

T ~ o o 0 x; jqy, yf 

0 0 0 2XjX 2 Xjy 2 +X 2 y, 2 y,y 2 

0 0 0 x\ x 2 y 2 y 2 2 _ 

One can also show that 

dV = dxdy = 2Vd£dri, ^ 3 ' 12 * 

where V = 1 1 r, x r, | is the volume of S. Substituting Eqs. (3. 1 1 ) and (3. 12) into Eq. (3.8), we 

obtain 
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'J c e^,n)d^dr\ ••• J c e N {£,n)dZdi) 

L = [e l (%,'n),...,e N (Z,ri)] : — '■ 

[ e^,r])d^di 1 ••• f e N &n)dtdr\ 

JC N JC .N 

From the above equation, the cardinal basis functions can be made independent of the nodal 
positions of S if each and every Vj is proportional to V. This can be achieved by subdividing the 
SV into polygonal CVs with straight edges. Therefore, different shapes of triangles have the 
identical cardinal bases £,,(£, 7?) in the transformed space D if they are similarly partitioned into 

polygons. Transforming back to the physical space S, although the functions Lj(x,y) may be 
different for different triangles, their functional values L } (x rq ,y rq ) at similar points (points 
having the same (£,fj) in the transformed space D) are exactly the same. We thus have a 
universal reconstruction formula, Eq. (3.9), for evaluating the state variable u at similar points. 
This also implies that the reconstruction needs to be carried out only once, and that can be 
performed using any shape of triangle. Although matrix R may be ill-conditioned, we avoid 
numerically inverting the matrix by using Mathematica [39] to derive the reconstruction 
coefficients analytically using exact arithmetic. These coefficients are identical for all triangles. 
The exact integrations of polynomials over arbitrary polygons can be found in [28]. 

Note that one of the subtle differences between a FV method and a SV method is that all the CVs 
in a SV use the same data reconstruction. As a result, it is not necessary to use a Riemann flux or 
flux splitting for the interior boundaries between the CVs inside a particular SV because the state 
variable is continuous across the interior CV boundaries. Riemann fluxes are only necessary at 
the boundaries of the SV. The most significant advantage of the SV method, as compared with the 
FV method, is that the reconstruction for a particular cell type (e.g. triangles) with a certain CV 
subdivision (e.g. those shown in Figures 2-4) is exactly the same. Therefore, the memory and 
CPU intensive reconstructions used in a FV method are solved analytically without taking any 
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extra memory in the SV method. Furthermore, exact fluxes rather than Riemann fluxes are used 
at the interior boundaries of the CVs, resulting again significant savings because the Riemann 
flux is usually several time more expensive to compute than the exact flux. 

4. Convergent Linear, Quadratic and Cubic Triangular Spectral Volumes 

Based on the discussions in the last section, it is clear that the reconstruction problem is 
equivalent for all triangles. We therefore focus our attention on the reconstruction problem in an 
equilateral triangle E, as shown in Figure 5b. In partitioning E into N non-overlapping CVs, we 
further require that the CVs satisfy the following three conditions: 

1. The CVs are "symmetric" with respect to all symmetries of the triangle; 

2. All CVs are convex; 

3. All CVs have straight sides, i.e., the CVs are polygons. 

We believe the symmetry and convexity requirement is important for achieving the best possible 
accuracy and robustness. The requirement of polygons simplifies the formulation of the SV 
method. To handle curved boundaries, isoparametric SVs will be used. The curved boundary will 
be represented using high-order polynomials compatible with the polynomial interpolation inside 
the SV. Then the isoparametric SVs will be transformed to the standard triangle, which will be 
partitioned in the usual manner. Surface integrals will be performed with respect to the standard 
triangle. It is obvious that a CV containing the centroid of E must be symmetric with respect to 
the three edges and vertices, and at most one such CV can exist. This CV, if it exists, is thus said 
to possess degree 1 symmetry (or 1 symmetry, in short). Similarly, CVs with degree 3 and 6 
symmetries can also be defined. For example, if a CV is said to possess degree 3 symmetry, then 
two other symmetric CVs must exist in the same partition. We shall denote n,, n 3 and n 6 the 
number of degree 1, 3 and 6 symmetry groups in a partition with n, = 0 or 1. Then the total 
number of CVs in the partition is then ny + 3n 3 + 6ti6. In order to support the unique 
reconstruction of a degree m polynomial, the total number of CVs must be identical to the 
dimension of the polynomial space, i.e.. 
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(4.1) 


(m + l)(m + 2) 

n, + 3n 3 + 6 n 6 = 

The solutions of (4.1) can be used to guide the partition of E once m is given. For example, for m 

- \ 5 > all possible solutions are summarized in Table 4.1. Some possible partitions of the 

standard triangle corresponding to these solutions for m = 1, 2, 3 are shown in Figures 2-4. Next 
the question of how these partitions perform in a data reconstruction needs to be answered. 

Given any partition, the reconstruction matrix R in (3.5) must be non-singular. In this case, the 
expansion coefficients can be solved from (3.6). Note that once the polynomial basis is given, the 
matrix is solely determined by the partition Tl m of E. For a linear reconstruction using three 

CVs, it is well known that the reconstruction is non-singular as long as the centroids of the CVs 
are not co-linear. Unfortunately no such simple criteria are known for higher order 
reconstructions in E. We will therefore have to compute the determinant of the reconstruction 
matrix to determine whether it is singular. As a matter of fact, straightforward computations 
indicated that several partitions shown in Figures 3 and 4 are singular. For example, the 
quadratic ST in Figure 3b and the cubic ST in Figure 4d are verified to be singular, and they will 
be excluded from further considerations. 

In [38], the first paper on the ST method, it was shown that not all non-singular reconstructions 
are convergent. For example, high-order polynomial reconstructions based on equidistant CVs in 
one dimension are not convergent although the reconstructions are non-singular. We believe this 
is the direct consequence of the Runge phenomenon. Therefore some means to quantify the 
quality of the reconstructions needs to be identified. 

Assume we have a non-singular partition FI m of E. For we C{E), we then have 

H 

The cardinal basis function has the following property 

— f LA€,Vi)d£dTi = 8 iJ , l<i,j<N. 

V Jc .- 


(4.2) 


(4.3) 
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Denote p m = r n («), where T n is an operator which maps C(E) onto P m (E) . It is obvious that 

r n is a linear projection operator because: 

• Tjj (m + v) = r n ( m) + Tpj (v), Vh 6 C{E), v & C(E) 

• r n (cM) = cr n (m) for any real constant c; 

• r n p = P for pe P m . 

When both spaces C{E) and PJE) are equipped with the supremum or uniform norm, i.e., 


l*ll = max|*| , the norm of this projection operator can be defined as 

IMhsupM 


u* 0 


( 4 . 4 ) 


Because u j — u \*j = , then it is obvious that 


ll^n w || _ 

7=1 

- < 

<- 

IHI 

7=1 

u 

U 


a 

HI 


Therefore we can easily see that 



= max 
§€£ 


7 = 1 


7=1 


( 4 . 5 ) 


( 4 . 6 ) 


The function A(£,t7) = £|L,(£,r?)| is usuall y referred to as the Lebesgue function of the 

>=l 

interpolation, and ||r n | is called the Lebesgue constant, which is of interest for the following two 
reasons [8]: 

(i) If P* m is the best uniform approximation to u on E, then 

t-r„«U(i+|Tn 1^-/4 


( 4 . 7 ) 


because 


u-T n u = 


\u- P m ~(Pn U ~ Prn )j| _ || M Pm Hi ( M Pm 


< 


« -Pi + r i 


— || M Pm + ll^n Pm ) | 

(ii) For peP m , 

< llr n || max u j. 


:«-p;)|=6+rni) 


“ - Pn 


1< j<N 


( 4 . 8 ) 
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which is obvious from (4.2). Thus ||r n || gives a simple method of bounding the interpolation 

polynomial. It is obvious from (4.7) that the smaller the Lebesgue constant, the better the 
interpolation polynomial is to be expected. Therefore the problem becomes finding the partition 
with small Lebesgue constant, if not as small as possible. In this paper, our focus is to construct 
good enough SV partitions so that the interpolation polynomial is convergent when the 
computational grid is refined. The Lebesgue constant is used as the criterion to judge the quality 
of the partitions. The optimization of the partitions will be the subject of a future publication. 

The problem of partitioning the equilateral triangle E into N symmetric CVs, which can support 
non-singular polynomial interpolations is not trivial, and is much more complex than the 
problem of determining a set of points in E which support Lagrange interpolations [8,1 1], What 
we try to accomplish in the paper is to identify partitions of E which support linear to cubic 
reconstructions with relatively small Lebesgue constants. 

Linear Spectral Volume {m - 1) 

Based on the solution of (4.1), it is obvious that two partitions are possible, as shown in Figure 
2a and 2b, which are named Type 1 and Type 2 partitions. Since the centroids of the CVs are 
non-co-linear, both partitions are admissible. Note that the CVs in both partitions possess a 
degree 3 symmetry. The cardinal basis functions L,(§,n) are plotted in Figure 6 for both the 

Type 1 and Type 2 partitions, respectively. Furthermore, the Lebesgue constants are 13/3 
(4.3333) and 43/15 (2.8667) for Type 1 and 2 partitions, respectively. Note that the Type 2 SV 
has a much smaller Lebesgue constant than the Type 1 SV, indicating that the error with the 
Type 2 SV should be smaller than the error with the Type 1 SV. 

For a linear reconstruction, only one Gauss quadrature point is required for a surface integral. 
This quadrature point is located at the center of an edge. Due to the symmetry, we therefore only 
need to compute and store the functional values of the cardinal bases at two quadrature points, 
total of six coefficients. These coefficients are the same for all triangles with similar partitions. 
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Quadratic Spectral Volume jm = 2) 

Eq. (4.1) has two solutions for m = 2, and the possible partitions are shown in Figure 3. As 
mentioned earlier, the partition shown in Figure 3b is singular, and is therefore not admissible. 
The partition presented in Figure 3a is not unique in the sense that the position of one of the two 
vertices on an edge of the triangle can change, i.e., the length d shown in Figure 3a can be any 
real number in (0, 0.5) assuming the length of the edge is 1. It seems that with any d, the partition 
is admissible. In our numerical studies, two different values of d were tested, namely d = 1/3 and 
d = 1/4 (corresponding to the Gauss-Lobatto points on the edge), which are called Type 1 and 
Type 2 partitions, respectively. The Lebesgue constant for the Type 1 partition is 9.3333, and for 
the Type 2 partition is 8. Therefore, the Type 2 partition is expected to yield more accurate 
numerical results. The cardinal basis functions L,(£,r?) are plotted in Figure 7 for the Type 2 

partition. For a quadratic reconstruction, two Gauss quadrature points are required for a surface 
integral. Due to the symmetry, there are 30 coefficients, corresponding the functional values of 
the cardinal bases at five quadrature points, which need to be computed and stored. 

Cubic Spectral Volume jm = 3} 

Eq. (4.1) has two different solutions for m = 3, and several possible partitions are shown in 
Figure 4. As mentioned earlier, the partition shown in Figure 4d is singular, and is not 
admissible. Although the partitions presented in Figure 4a-4c look quite different, the first two 
partitions can be viewed as the limiting cases of the partition shown in Figure 4c. Therefore we 
can claim that all three partitions 4a-4c have the same general topology. The partition requires 
the locations of three vertices I, II and III as shown in Figure 4c. The optimization of these points 
will be the topic of a future paper. In this paper, our focus is the partition shown in Figure 4b, in 
which one parameter d can be changed to obtain different partitions. In fact, the Lebesgue 
constants for partitions with a set of d values are presented in Table 4.2. Among this set of d 
values, it is interesting to note that the Lebesgue constant reaches a smallest value of 3.44485 at 
j - 1/15 from a value of 8.21499 at d = 1/6. When d is smaller than 1/15, the Lebesgue constant 
starts to increase. For presentation purpose, we call the partition shown in Figure 4a the Type 1 
partition. The partition shown in Figure 4b with d = 1/6 is called the Type 2 partition, and with d 
= 1/15 the Type 3 partition. It is expected that the Type 3 partition should give the most accurate 
numerical solution in the uniform norm. The Lebesgue constant for the Type 1 partition is 
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167/12 (13.9), which is significantly larger than those for the Type 2 and 3 partitions. Numerical 
results to be presented later confirm that the Type 1 partition is not convergent with grid 
refinement. Several of the cardinal basis functions for the Type 1 partition are plotted in Figure 
8, and for the Type 3 partition are plotted in Figure 9. For a cubic reconstruction, still two Gauss 
quadrature points are required for a surface integral. Again, due to the symmetry, we only need 
to compute and store the cardinal basis functional values at 10 quadrature points, total of 100 

coefficients. 


5. Multi-Dimensional Limiters 

The Gibbs phenomenon associated with high-order schemes in the presence of discontinuities 
causes loss of monotonicity in the solution of hyperbolic conservation laws. Godunov [18] first 
proved that there are no linear second or higher order schemes which guarantee monotonicity. 
Therefore high-order monotonic schemes, if they exist, must be non-linear. One effective 
approach to achieve monotonicity is to limit the reconstructed solution so that a monotonicity 
constraint is satisfied. As pointed out by many researchers, strict monotonicity seems to conflict 
with uniform high-order of accuracy [12]. In order to recover uniform order of accuracy away 
from discontinuities, the TVB (total variation bounded) idea [33] is employed here. 


We again consider a SV 5, with N CVs. Given cell-averaged state variables for all the CVs {u tJ } 


a polynomial reconstruction p,{x,y) of at most order k - 1 exists which satisfies. 

l Pi (x,y)dV = j = 


(5.1) 


A, 


Recall that this polynomial reconstruction is then used to compute the state variables at the CV 
boundaries, which are, in turn, used in the update of the solution at the next time level: 


K J 


7 — + ^ X ^Rkm ( P, ( X rq ’ )% )’ /Ar ( X rq ’ ? r,, )’ H r ) A- 

Clt V,,j r=l <7=1 


(5.2) 


Denote 
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A u rq = Pi (x rq , y rq ) ~ u Lj , r - 1,- • • , K\ q - 1, • ■ * , J- 
Following the TVB idea, if 

\tXu rq \<4Mh; q , r = l = ( 5 - 3 ) 

it is not necessary to do any data limiting. In (5.3), M represents some measure of the second 
derivative of the solution, and h rq is the distance from point (x rq , y rq ) to the centroid of C u . It is 
obvious that this multi-dimensional TVB limiter degenerates into the one-dimensional TVB 
limiter. Using the fact that in two dimensions, 4h 2 q V i } , we can further simplify (5.3) by 

replacing 4Mh 2 q with MV uj • With the new formula, we do not need to compute or store the 

distances from the cell centroid to the quadrature points. In this paper, we select M to be close to 
the maximum absolute value of the second derivative over the computational domain. If for any 
value of r and q, (5.3) is violated, it is assumed that C u is near a steep gradient and data limiting 
is necessary. Instead of using the polynomial Pi (x, y ) in C u , we assume that data is linear m Cy, 


i.e. 


u l j (x,y) = u l j +Vm t j •(r-r ij ), Vre C t J , 


(5.4) 


where r t ; is the position vector of the centroid of C u . In order to achieve the highest resolution, 
we need to maximize the magnitude of the solution gradient Vm,- j in C,j. At the same time, we 
require that the reconstructed solutions at the quadrature points of C u satisfy the following 
monotonicity constraint: 




<u,,j(x rq ,y rq )<u^\ r 


(5.5) 


where n/™ 1 and u™-* are the minimum and maximum cell-averaged solutions among all its 

neighboring CVs sharing a face with Cy, i.e., 

u m r = max(M, . , max m, , r ) 

‘■J '■> 1 < r <K J 


U m0 

U ‘.J 


= min (U: ,, min u, , ), 

l<r<X 1 


(5.6) 
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where u Lj _ r denotes the cell-averaged solution at the neighboring CV of Cy shanng face r. If the 
solution is linear, it is obvious that the solution value at the centroid of C u is the same as u iJt 
the cell-average solution. Several different approaches are possible in estimating Vn,j 
depending on the computational complexity. Three approaches are outlined here. 


Approach 1: 

Consider a CV with K faces. Using cell-averaged data at Cy and its neighbors, we can usually 
construct K different gradients. For example, the gradient reconstructed from 
1) ls denoted by V«„- r with 5f itM+i = u iJA • In addition, another gradient 

can be constructed through a least squares reconstruction algorithm using the cell-averaged data 
at all the face-neighbor cells. This gradient is denoted by Vu iJK+l . Using any of the gradients, 
the state variable at the quadrature points of Cy can be computed. If any of the reconstructed 
variable at the quadrature points is out of the range [ u L j , w, j ]. ih e gradient is limited, i.e., 

Vn„ ; , r «=<pV M(Jr , (5 ' 7) 

where (p £ [0, 1] is calculated from. 


V 


miru 


1 , 


An 


rq 


—max — 
u i,j u i,j 


if A U rq > 0 


min] 

1 


__ All rq 


V 


— —min 


if A u rq < 0. 
otherwise 


(5.8) 


After this limiting step, we have K+J gradients which all satisfy the monotonicity constraint 
given in (5.5). Then the gradient with the largest magnitude is selected to maximize the gradient, 


i.e., 


Vu, j = max Vn ( y r 


(5.9) 
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This limiter is similar to the Maximum Limited Gradient (MLG) limiter developed by Batten, 
Lambert and Causon [7], It was proven by Liu [26] that the above limiter satisfies a maximum 
principle for triangular grids. The MLG limiter is an analog of the Superbee limiter in one 
dimension. This limiter is therefore called Superbee limiter in this paper. This limiter has the 
advantage of minimum numerical dissipation, but is expensive to compute. A more efficient 

limiter is given next. 

Approach 2: 

Only one gradient is computed, i.e., V M ,„* +1 given in Approach 1, which is computed with data 

at all face-neighbor cells using a least squares linear reconstruction algorithm. Again this 
gradient is limited using (5.8) so that the reconstructed solutions at all the quadrature points 
satisfy the monotonicity constraint (5.5). This limiter is similar to the minmod limiter in one 
dimension, and is also called the Minmod limiter here. 

Approach 3: 

Note that in both approaches 1 and 2, the gradients of the solution for each CV must be 
reconstructed using data in a neighborhood of the CV. This reconstruction can be quite memory 
and CPU intensive. In this approach, we avoid a separate data reconstruction by reusing the 
polynomial reconstruction already available for the SV. For each CV, we use the gradient of the 
reconstructed polynomial at the CV centroid, i.e., 

_ Jdp m + + . (5.10) 

u,j dx ’ dy l [ dt, dx dr] dx’ dy drj dy 

For the transformation given in (3.10), it is obvious that 

dp m } f d Pm . ( OEjC 

17 = 1 _ * -*) * . (5.11) 

fynt 7 V y 1 dPm -V,y 2 - X 2)'li _X 2 X 1 J 

^ J I ^ j ^ 
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Because 


7=1 

we then have 

dpj^/n) _ V (5.12a) 

K 11 

^pj^n} = y ( 5 . 12 b) 

at? £ an 

The first derivatives of the shape function can be obtained analytically. The gradient for each CV 
is then limited if necessary with the same approach outlined in Approach 1. Obviously this is 
most efficient among the three limiters. This limiter is named CV limiter. 

For comparison purposes, we also used another simple limiter, which is called the Clip limiter. 
In this limiter, zero gradient (piece-wise constant distribution) is used in CVs wherever limiting 

is necessary. 

Note that if parameter M = 0, the TVB limiters are similar to TVD (total variation diminishing) 
limiters, which strictly enforce monotonicity by sacrificing accuracy near extrema. 

The availability of cell-averaged data on the CVs inside a SV makes this CV-based data limiting 
possible, whereas in the DG method, one can only do an element based data limiting. Due to the 
increased local resolution, the SV method is expected to have higher resolutions for 
discontinuities than the DG method. The improved resolution has been demonstrated in one 

dimension [38]. 
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6. Numerical Tests 


1. Accuracy Study with 2D Linear Wave Equation 

Time Accurate Problem 

In this case, we test the accuracy of the SV method on the two-dimensional linear equation: 

-1<*<1, - 1 < y < L 61) 

dt dx 3y ' ‘ ’ 

u (x, y, 0) = u 0 (x, y ), periodic boundary condition. 

The initial condition is u 0 (*,y) = sm*(* + y) ■ A fourth-order accurate Gauss quadrature 

formula [28] is used to compute the CV-averaged initial solutions. These CV-averaged solutions 

are then updated at each time step using the third-order TVD Runge-Kutta scheme presented 

earlier. The numencal simulation is earned until t = 1 on two different tnangular gnds as shown 

in Figure 10. One gird is generated from a uniform Cartesian grid by cutting each Cartesian cell 

into two triangles, and is named the regular grid. The other gnd is generated with an unstructured 

grid generator, and is named the irregular grid. The finer irregular grids are generated recursively 

by cutting each coarser grid cell into four finer grid cells. Note that the cells m the irregular grid 

have quite different sizes. In Table 6.1, we present the L, and errors in the CV-averaged 

solutions produced using second to fourth order SV method schemes with SVs shown in Figures 

2-4 on the regular gnd. The errors presented in the table are time-step independent because the 

time step At was made small enough so that the errors are dominated by the spatial 

discretization. In this test, all SVs except the Type 1 cubic SV (shown in Figure 4a) are 

convergent with grid refinement on this regular gnd. It is obvious that the expected order of 

accuracy is achieved by all the convergent SVs. It is not surprising that the Type 1 cubic SV is not 

convergent because of its rather large Lebesgue constant of 13.9. In contrast, the Type 2 and 3 

cubic SVs have Lebesgue constants of 8.21 and 3.44 respectively. It is interesting to note that the 

Type 1 linear SV gives more accurate results in both the L, and norms than the Type 2 linear 

SV even if the Type 1 S V has a larger Lebesgue constant of 4.33 versus that of 2.87, of the Type 

2 SV. This indicates that the Lebesgue constant cannot serve as an absolute error estimate, but 

rather an estimate of the upper bound of the error. For the quadratic and cubic SVs, the partitions 

with smaller Lebesgue constants do give more accurate numerical solutions, as shown in Table 

6 . 1 . 


23 


Next, the L[ and L errors in the numerical results computed on the irregular grid using second 
fourth order SVs are shown in Table 6.2. This should be a much tougher test case because of the 
truly unstructured nature of the computational grid. What is striking is that the Type 1 linear S\ 
failed to achieve second-order accuracy on this grid. As a matter of fact, it is only first order 
accurate. This may be contributed to the acute angles of the CVs in the partition. Note that both 
quadratic SVs are convergent, and give similar results. Third order accuracy is achieved by both 
types of quadratic SVs in the L, norm although the numerical order of accuracy in the L norm 
is only slightly over second-order. We believe this is due to the non-smoothness of the 
computational grid. As expected, the Type 1 cubic SV is not convergent on this grid. In addition, 
the Type 2 cubic SV also showed a non-convergent behavior in the norm on the finest grid. It 
is nice to see that the Type 3 cubic SV is not only convergent, but also fourth-order accurate in 
both the L, and norms. 

Steady State Problem 

One steady state solution for the wave equation (6.1) is u(x,y) = sin7r(x- y) . In order to test the 
performance of the SV method for steady-state problems, the steady boundary value problem is 
also studied. Because the wave is traveling in positive x and y directions, inflow boundary 
conditions are employed at y = -1 and x = -1, while extrapolation boundary conditions are used 
at x = 1 and y = 1. On the inflow boundaries, the exact solutions at the quadrature points are 
used to evaluate the flux integrals. The simulation is carried out on the irregular grid until a 
steady-state is reached. In all the simulations, the residuals were reduced to machine zero. In 
Table 6.3, the L[ and L„ errors are presented for second to fourth-order schemes. It is surprising 
to see that the Type 1 linear SV gives a more accurate solution in both norms than the Type 2 
linear SV on this irregular grid. Recall that in the time accurate simulation presented earlier on 
this irregular grid, the Type 1 linear SV failed to achieve second-order accuracy. This may 
indicate that there is significant error accumulation in the time accurate simulation with the Type 
1 linear SV. Other than that, there are no major surprises. Both quadratic SVs gave reasonable 
results, while the Type 1 and Type 2 cubic SVs showed convergence problems on the finest 
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mesh. Once again, the performance of the Type 3 cubic SV is excellent. It is convergent, and 
achieves fourth-order accuracy in both norms. 


2. Accuracy Study with 2D Burgers Equation 

In this case, we test the accuracy of the SV method on the two-dimensional non-linear wave 
equation: 

du 3m 2 12 ( dtr_[2 _!<*<!, -l<y<l, 

dt dx dy (6.2) 

n (x v 0) = - + - sin n{x + y), periodic boundary condition. 

4 2 

The initial solution is smooth. Due to the non-lineanty of the Burgers equation, discontinuities 
will develop in the solution. Therefore we also test the capability of the SV method to achieve 
uniform high-order accuracy away from discontinuities. At t = 0.1, the exact solution is still 
smooth, as shown in Figure 11a. The numerical simulation is therefore earned out until t = 0.1 
without the use of limiters on both the regular and irregular grids as shown in Figure 10. The 
numerical solution on the 20x20x2 irregular gnd computed with the Type 1 quadratic SV (third- 
order accurate) is displayed in Figure lib. Notice that the visual agreement between the 
numerical and exact solutions is excellent. In Table 6.4, we present the L, and L x errors 
produced using vanous SVs on the regular gnd, while in Table 6.5 the errors on the inegular 
grids are presented. The Type 1 cubic SV is now excluded because it is non-convergent on any 
grids. The performance of the SV method on the non-linear Burgers equation is quite similar to 
the performance on the linear wave equation, although there is a slight loss of accuracy (from 0. 1 
- 0.6 orders) especially on the inegular grid in the norm, probably due to the non-linear 
nature of the Burgers equation. Once again, the Type 1 linear SV has difficulty in achieving 
second-order accuracy on the inegular grid in both norms. 


At t = 0.45, the exact solution has developed two shock waves as shown in Figure 12a. Limiters 
are necessary to handle the discontinuities. All the limiters (Clip, CV, Superbee and Minmod 
limiters) are evaluated. Shown in Figure 12 are the exact solution, and the computed numerical 
solutions with the Type 2 quadratic SV (third order accurate) on the 40x40x2 inegular gnd using 
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all limiters with M = 0, i.e., TVD limiters. The use of M - 0 was designed to highlight the 
differences between the limiters. Note that all the TVD limiters gave reasonable solutions 
including the Clip limiter. In term of shock resolution, the Clip Limiter is the most dissipative, 
followed by the CV Limiter. The Minmod and Superbee Limiters gave the best solutions, which 
are difficult to distinguish from each other. Figure 13 displays the solutions with the same grid 
and SV scheme with M = 400, i.e., TVB limiters. The most striking difference between the 
results shown in Figure 12 and 13 is that the solutions away from the shock waves are much 
smoother with the TVB limiters than with the TVD limiters. This indicates that high-order 
accuracy is achieved away from the discontinuities with the TVB limiters. However, the price 
one must pay to achieve this is that some oscillations near the shock waves must be tolerated, as 

shown in Figure 13. 

In order to estimate the numerical order of accuracy for the solution away from the 
discontinuities, L, and L x errors in the smooth region [-0.2, 0.4]x[-0.2, 0.4] are computed. 
Computations were carried out on the irregular grid only with the Minmod limiter. Without the 
use of the limiter, the solution quickly diverged after shock waves were developed in the 
solution. The parameter M was set to be 400 in the computation. If M is too small, the accuracy 
in the smooth region is degraded probably because limiting was carried out in the smooth region 
as well as near the shock. The L, and errors with one type of SV for a given order of 
accuracy are presented in Table 6.6. Obviously, with this choice of M, the designed order of 
accuracy was achieved away from discontinuities. 

7. Conclusions 

The Spectral Volume method [38] has been successfully extended to two-dimensional scalar 
conservation laws using unstructured triangular meshes. Each mesh cell forms a spectral volume, 
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and the spectral volume is further partitioned into polygonal control volumes. High order 
schemes are then built based on the CV-averaged solutions. It is shown that a universal 
reconstruction can be obtained if all spectral volumes are partitioned in a similar manner. 
However, as in the one-dimensional case, the way in which a SV is partitioned into CVs affects 
the convergence property of the resultant numerical scheme. A criterion based on the Lebesgue 
constant has been developed and used successfully to determine the quality of various partitions. 
Symmetric, stable, and convergent linear, quadratic and cubic SVs have been obtained, and many 
different types of partitions are evaluated based on the Lebesgue constants and their performance 

on model test problems. 

Accuracy studies with 2D linear and non-linear scalar conservation laws have been earned out, 
and the order of accuracy claim has been numerically verified on both smooth and non-smooth 
triangular gnds for convergent SVs. Several TVD and TVB limiters have been developed for 
non-oscillatory captunng of discontinuities, and found to perform well. The TVB limiters with a 
properly selected parameter (A#) are capable of maintaining uniformly high-order accuracy away 
from discontinuities. The extension of the method to one and two dimensional hyperbolic 
systems is under way, and will be reported in future publications. 
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Table 4.1 Solutions of (4.1) for m = 1, ...,5 
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Table 6.1. Accuracy on u t + u x + u y 


= 0, with uq (jc, y ) = sin n(x + y), at t = 1 (regular grids) 


Order of Accuracy 

Grid 

L] error 

L[ order 

L,„ error 

L, x order 


10x10x2 

3.04e-2 

- 

4.97e-2 

- 


20x20x2 

7.68e-3 

1.98 

1.24e-2 

2.00 

2 

40x40x2 

1.92e-3 

2.00 

3.10e-3 | 

2.00 

(Type 1 5 VO 

80x80x2 

4.81e-4 

2.00 | 

7.75e-4 

2.00 

160x160x2 

1.20e-4 

2.00 

1.93e-4 

2.00 ! 

2 

(Type 2 SV) 

10x10x2 

4.03e-2 

- 

6.68e-2 


20x20x2 

1.06e-2 

1.93 

1.78e-2 

1.91 

40x40x2 

2.7 le-3 1 

1.97 

4.54e-3 

1.97 1 

80x80x2 

6.83e-4 

1.99 

1.14e-3 » 

1.99 

160x160x2 

1.71e-4 

2.00 

2.87e-4 

1.99 

3 

(Type 1, d = 1/3) 

10x10x2 

4. 18e-3 

- 

7.76e-3 

- 

20x20x2 

5.33e-4 

2.97 

1.01e-3 

2.94 

40x40x2 

6.73e-5 

2.99 

1.25e-4 

3.01 

80x80x2 

8.45e-6 

2.99 

1.55e-5 

3.01 

160x160x2 

1.06e-6 

2.99 

1.93e-6 

3.00 

3 

(Type 2, d = 1/4) 

10x10x2 

4.73e-3 

- 

7.88e-3 

- 

20x20x2 

4.77e-4 

3.31 

9.83e-4 

3.00 

40x40x2 

6.04e-5 

2.98 

1.23e-4 

3.00 

80x80x2 

7.58e-6 

2.99 

1.53e-5 

3.01 

160x160x2 

9.57e-7 

2.99 1 

1.91e-6 

3.00 

4 

(Type 1 SVO 

10x10x2 

1.38e-4 


4.86e-4 

- 

20x20x2 

8.64e-6 

4.00 

1.98e-5 

4.62 

40x40x2 

5 Ale-7 

3.98 

1.51e-6 

3.71 

80x80x2 

3.46e-8 

3.98 

1.17e-7 

3.69 

160x160x2 

4.19e-8 

Negative 

5.15e-7 

Negative 

4 

(Type 2, d = 1/6) 

10x10x2 

9.33e-5 

- 

3.17e-4 

" 

20x20x2 

5.86e-6 

3.99 

1.94e-5 

4.03 

40x40x2 

3.70e-7 

3.99 

1.24e-6 

3.95 

80x80x2 

2.32e-8 

4.00 

7.78e-8 

3.99 

160x160x2 

1.45e-9 

4.00 

4.84e-9 

4.01 

4 

(Type 3, d = 1/15) 

10x10x2 

7.36e-5 

- 

2.51e-4 

" 

20x20x2 

4.52e-6 

4.03 

1.61e-5 

3.96 

40x40x2 

2.81e-7 

4.01 

1. Ole-6 

3.99 

80x80x2 

1.75e-8 

4.01 

6.30e-8 

4.00 

160x160x2 

1.10e-9 

3.99 

3.94e-9 

4.01 
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Table 6.2. Accuracy on u, + u x + u y 


= 0, with u 0 (x, y ) = sin^(x + y), at t — l (irregular grids) 


Order of Accuracy 

Grid 

L} error 

Lj order 

L,„ error 

order 


10x10x2 

1.30e-l 1 

- 

3.60e-l 

- 

2 

20x20x2 

6.66e-2 

0.96 

1.91e-l 

0.91 

(Type 1 SV) 

40x40x2 

3.51e-2 

0.92 

9.84e-2 [ 

0.96 


80x80x2 

1.85e-2 

0.92 ~T 

4.91e-2 

1.00 1 


160x160x2 

9.74e-3 

0.93 

2.86e-2 

0.78 


10x10x2 

6.71e-2 

- 

1.36e-l 

- 

2 

(Type 2 SV) 

20x20x2 

1.83e-2 

1.87 

4.42e-2 f 

1.62 

40x40x2 

4.71e-3 

1.96 

1.15e-2 

1.94 

80x80x2 

1.19e-3 

1.98 

2.94e-3 

1.97 

160x160x2 

3.00e-4 

1.99 

8.85e-4 

1.73 

3 

(Type 1, d = 1/3) 

10x10x2 

9 Ale-3 

- 

3.67e-2 

- 

20x20x2 

1.25e-3 

2.87 

5.28e-3 

2.80 

40x40x2 

1 ,64e-4 

2.93 

8.32e-4 

2.67 

80x80x2 

2.15e-5 

2.93 

1.84e-4 

2.18 

160x160x2 

2.79e-6 

2.95 

4.05e-5 

2.18 

3 

(Type 2, d = 1/4) 

10x10x2 

8.36e-3 

- 

3.76e-2 


20x20x2 

1.15e-3 

2.86 

5.63e-3 

2.74 

40x40x2 

1.52e-4 

2.92 

1.00e-3 

2.49 

80x80x2 

2. 01e-5 

2.92 

2.14e-4 

2.22 

160x180x2 

2.64e-6 

2.93 

5.31e-5 

2.01 

4 

(Type 1 SV) 

10x10x2 

4.43e-4 

- 

3.32e-3 

“ 

20x20x2 

3.08e-5 

3.85 

2.56e-4 

3.70 

40x40x2 

2.15e-6 

3.84 

2.12e-5 

3.59 

80x80x2 

2.48e-7 

3.12 

5.06e-6 

2.07 

160x160x2 

5.19e-7 

Negative 

4.49e-5 

Negative 

4 

(Type 2, d = 1/6) 

10x10x2 

3.04e-4 

- 

2.58e-3 


20x20x2 

2.02e-5 

3.91 

1.73e-4 

3.90 

40x40x2 

1.34e-6 

3.91 

1.42e-5 

3.61 

80x80x2 

9.61e-8 

3.80 

1.03e-6 

3.79 

160x160x2 

2.30e-8 

2.06 

1.23e-6 

Negative 

4 

(Type 3, d= 1/15) 

10x10x2 

2.71e-4 

- 

1.51e-3 

- 

20x20x2 

1 .6 le-5 

4.07 

1.14e-4 

3.73 

40x40x2 

9.91e-7 

4.02 

8.28e-6 

3.78 

80x80x2 

6.17e-8 

4.01 

5.40e-7 

3.94 

160x160x2 

3.87e-9 

3.99 

3.79e-8 

3.83 
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Table 6.3. Accuracy on u, + u, = 0, with u(x. y) = sin rr(.r - y) (irregular gods) 


Order of Accuracy 

Grid 

L\ error 


L ot error 



10x10x2 

8.91e-3 

- 

4.65e-2 

- 

2 

20x20x2 

2.26e-3 

1.98 

1.27e-2 

1.87 

(Type 1 SV) 

40x40x2 

5.78e-4 

1.97 

3.37e-3 

1.91 


80x80x2 

1.48e-4 

1.97 

9.03e-4 

1.90 

" 

160x160x2 

3.77e-5 

1.97 

2.37e-4 

1.93 


10x10x2 

1.57e-2 

- 

6.58e-2 

- 

2 

20x20x2 

3.94e-3 

1.99 

2.14e-2 

1.62 

(Type 2 SV) 

40x40x2 

9.99e-4 

1.99 

5.58e-3 

1.94 


80x80x2 

2.52e-4 

1.99 

1.47e-3 

1.92 j 


160x160x2 

6.32e-5 

2.00 

3.83-4 

1.94 1 

3 

(Type 1, d= 1/3) 

10x10x2 

1 .5 le-3 

- 

1.33e-2 

- 

20x20x2 

2.20e-4 

2.78 

1.82e-3 

2.87 ' 

40x40x2 

3.07e-5 

2.84 

2.64e-4 

2.75 

80x80x2 

4.08e-6 

2.91 

3.93e-5 

2.75 

160x160x2 

5.24e-7 

2.96 

5.25e-6 

2.90 

3 

(Type 2, d = 1/4) 

10x10x2 

1.48e-3 

- 

1.26e-2 

“ 

20x20x2 

2.1 le-4 

2.81 

1.72e-3 

2.87 

40x40x2 

2.88e-5 

2.87 

2.7 le-4 

2.67 

80x80x2 

3.76e-6 

2.94 

3.83e-5 

2.82 

160x160x2 

4.69e-7 

3.00 

5.11e-6 

2.91 

4 

(Type 1 SV) 

10x10x2 

9.11e-5 

- 

1.28e-3 

“ 

20x20x2 

6.80e-6 

3.74 

1 .04e-4 

3.62 

40x40x2 

4.21e-7 

4.01 

7.60e-6 

3.77 

80x80x2 

2.97e-8 

3.83 

8.75e-7 

3.12 

160x160x2 

3.73e-9 

2.99 

7.21e-7 

0.28 

4 

(Type 2, d = 1/6) 

10x10x2 

9.27e-5 

- 

1.26e-3 

“ 

20x20x2 

6.34e-6 

3.87 

1.02e-4 

3.63 

40x40x2 

3.93e-7 

4.01 

7.36e-6 

3.79 

80x80x2 

2.56e-8 

3.94 

4.83e-7 

3.93 

160x160x2 

1.76e-9 

3.86 

6.25e-8 

2.95 

4 

(Type 3, d = 1/15) 

10x10x2 

1.59e-4 

- 

9.2 le-4 

- 

20x20x2 

9.98e-6 

3.99 

9.12e-5 

3.34 

40x40x2 

6.26e-7 

3.99 

6.58e-6 

3.79 

80x80x2 

3.75e-8 

4.06 

4.49e-7 

3.87 

160x160x2 

2.26e-9 

4.05 

2.68e-8 

j[ 4.07 
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Table 6.4. Accuracy on u, + uu x + uu y - 0, with u 0 (x, y) - ^ + 

regular grid 


— sin7r(x+ y ), at t 
2 


Order of Accuracy 

Grid 

error 


10x10x2 

2.05e-3 

9 

20x20x2 

5.31e-4 


40x40x2 

1.42e-4 

(Type 1 SV) 

80x80x2 

3.79e-5 


160x160x2 

1.03e-5 


10x10x2 

3.35e-3 


20x20x2 

8.07e-4 


40x40x2 

2.01e-4 

(Type 2 SV) 

80x80x2 

5.06e-5 


160x160x2 

1.27e-5 


10x10x2 

4.44e-4 


20x20x2 

7.81e-5 

3 

(Type 1, d = 1/3) 

40x40x2 

1.30e-5 

80x80x2 

2.09e-6 


160x160x2 

3.24e-7 


10x10x2 

4.31e-4 

3 

(Type 2, d = 1/4) 

20x20x2 

7.44e-5 

40x40x2 

1.23e-5 

80x80x2 

1.96e-6 


160x160x2 

3.00e-7 


10x10x2 

3.26e-5 

A 

20x20x2 

2.22e-6 

4 

(Type 2, d = 1/6) 

40x40x2 

1.58e-7 

80x80x2 

1.02e-8 


160x160x2 

6.54e-10 


10x10x2 

4.10e-5 

A 

20x20x2 

2.69e-6 

4 

(Type 3, d = 1/15) 

40x40x2 

1.85e-7 

80x80x2 

1.24e-8 


160x160x2 

8.21e-10 


L] order 





L„ error 

7.20e-3~ 


2.17e-3 


6.63e-4 


1.96e-4 


5.86e-5 


1.18e-2 


3.82e-3 


1.02e-3 


2.64e-4 


6.71e-5 


1.93e-3 


4.29e-4 


8.84e-5 


1.51e-5 


2.57e-6 


1.91e-3 


4.22e-4 


7.93e-5 


1.45e-5 


2.42e-6 


2.40e-4 


2.15e-5 


1.74e-6 


1.18e-7 


7.68e-9 


3.12e-4 


2.82e-5 


2.23e-6 


1.60e-7 


1.09e-8 
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Table 6.5. Accuracy on u, + uu x + uu y - 0, with u 0 (x,y) - ^ 

with irregular grid 


+ — sin7r(jt + v), at t = 0.1 
2 


Order of Accuracy 


(Type 1 SV) 


(Type 2 SV) 


(Type 1, d = 1/3) 


(Type 2, d = 1/4) 


(Type 2, d = 1/6) 


(Type 3, d = 1/15) 


10x10x2 


20x20x2 


40x40x2 


80x80x2 


160x160x2 


10x10x2 


20x20x2 


40x40x2 


80x80x2 


160x160x2 


10x10x2 


20x20x2 


40x40x2 


80x80x2 


160x160x2 


10x10x2 


20x20x2 


40x40x2 


80x80x2 


160x160x2 


10x10x2 


20x20x2 


40x40x2 


80x80x2 


160x160x2 


10x10x2 


20x20x2 


40x40x2 


80x80x2 


160x160x2 I 2.88e-9 



5.02e-7 


6.28e-4 


1.17e-4 


1.91e-5 


3.01e-6 


4.63e-7 


7.87e-5 


6.07e-6 


4.55e-7 


3.44e-8 


2.79e-9 


9.71e-5 


7.17e-6 


5.20e-7 


3.79e-8 
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Table 6.6. Accuracy on u, + uu x + uu y = 0, with u 0 (x, y ) - - + -sin n(x + y), at t 0.45 
in [-0.2, 0.4]x[-0.2, 0.4] on irregular grid, Minmod Limiter with M = 400 


Order of Accuracy 

Grid 

L { error 

Ly order 

L„ error 

L„ order 

2 

(Type 2 SV) 

10x10x2 

1.68e-4 

- 

5.33e-3 

~ 

20x20x2 

3.92e-5 

2.10 

1.65e-3 

1.69 

40x40x2 

9.66e-6 

2.02 

4.83e-4 

1.77 

80x80x2 

2.43e-6 

1.99 

1.58e 4 

1.61 

160x160x2 

6. 01e-7 

2.02 

3.49e-5 

2.18 

3 

(Type 2 SV) 

10x10x2 

6.23e-5 

- 

6.57e-3 

“ 

20x20x2 

6.25e-6 

3.32 

5.86e-4 

3.49 

40x40x2 

6.06e-7 

3.37 

7.21e-5 

3.02 

80x80x2 

7.40e-8 

3.03 

1.29e-5 

2.48 

160x160x2 

9.47e-9 

2.97 

2.68e-6 

2.27 

4 

(Type 3 SV) 

10x10x2 

7.81e-5 

- 

4.39e-2 

- 

20x20x2 

6.78e-7 

6.85 

1.3 le-3 

5.07 

40x40x2 

6.38e-9 

6.73 

2.97e-6 

8.78 

80x80x2 

3.85e-10 

4.05 

6.65e-8 

5.48 

160x160x2 

2.84e-l 1 

3.76 

4.36e-9 

3.93 
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Figure 1. (a) A possible reconstruction stencil for a quadratic reconstruction in a high-order k- 
exact finite volume scheme; (b) The partition of a Spectral Volume into six Control Volumes 

supporting a quadratic reconstruction 
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(a) Type 1, rij = 0, 113 - 1, n6 - 0 



(b) Type 2, ni = 0, 113 - 1, n6 - 0 

Figure 2. Control Volumes in a Triangular Linear Spectral Volume (a) Type 1 ; (b) Type 2 




(a) Type 1 with d = 1/3 and Type 2 with d - 1/4 
m = 0, n 3 = 2, n 6 = 0 



(b) A Singular Partition 
ni = 0, n 3 = 0, n<, = 1 


Figure 3. Possible Triangular Quadratic Spectral Volume Partitions 
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(b) Type 2 with d = 1/6 and Type 3 with d - 1/15 
m = l,n 3 = l,n 6 = 1 




Figure 4. Possible Cubic Triangular Spectral Volumes 
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(a) 



(b) 

Figure 5. The Schematic of the mapping from the physical triangle to the standard mangle 
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(b) Type 2 

Figure 6. Shape Functions for Linear Spectral Volumes 
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(b) 

Figure 7. Shape Functions for the Type 2 Quadratic Spectral Volume 
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45 
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(a) Exact Solution at t = 0. 1 


(b) Numerical Solution on the 20x20x2 Irregular Grid Using the Type 1 Quadratic SV 
Figure 11. Exact and Computational Solutions of the Burgers Equation at t = 0. 1 
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(a) Exact Solution (b) Clip Limiter 



(e) Superbee Limiter 

'igure 12. Exact and Computational Solutions of the Burgers Equation at t = 0.45 on the 
40x40x2 Irregular Grid Using the Type 2 Quadratic SV (Third-Order Accurate), M - 0 









(c) CV Limiter (d) Minmod Limiter 



(e) Superbee Limiter 

Figure 13. Exact and Computational Solutions of the Burgers Equation at t = 0.45 on the 
40x40x2 Irregular Grid Using the Type 2 Quadratic SV (Third-Order Accurate), M = 400 








